function [pop_cf,lambda_cf,L_k_cf,L_y_cf,L_o_cf,pi_cf] = counter_execute(N,S,d_geo,d_know,d_geo_distance,d_migr_const_aux,migration_exposure_cf,intra_id,lambda_start,alpha_cf,epsilon_cf,...
    pop_start,L_k_start,L_y_start,ubar_cf,mu,theta,zeta,omega,fertility_cf)

%%% Define diffusion frictions:

diffusion_frictions_cf = zeros(N*S,N*S);

d_migr = ones(N,N);
for i=1:N
    for j=1:N
        d_migr(i,j) = exp(d_migr_const_aux*migration_exposure_cf(i,j));
    end
end

for i=1:N
    for s=1:S 
        for j=1:N
            for r=1:S
                if intra_id == 0
                    diffusion_frictions_cf((s-1)*N+i,(r-1)*N+j) = (1/(d_geo(i,j)*d_migr(i,j)*d_geo_distance(i,j)*d_know(s,r)))^theta;
                else
                    if r == s
                        diffusion_frictions_cf((s-1)*N+i,(r-1)*N+j) = (1/(d_geo(i,j)*d_migr(i,j)*d_geo_distance(i,j)))^theta;
                    else
                        diffusion_frictions_cf((s-1)*N+i,(r-1)*N+j) = 0;
                    end
                end
            end
        end
    end
end

%%% Predict lambdas:

Summation_cf = zeros(N,S);
for j=1:N
    for r=1:S
            Summation_cf(j,r) = sum(reshape(lambda_start(:,:),N*S,1).*...
                diffusion_frictions_cf(:,(r-1)*N+j).*reshape(repmat(alpha_cf'.^theta,N,1),N*S,1)); 
            % Note the transpose on alpha_cf'. This is because alpha_cf is
            % a Sx1 vector. We want to transpose it into a 1xS vector, then
            % copy it N times vertically, so to get a NxS matrix
    end
end

lambda_cf = lambda_start + (epsilon_cf.^theta).*Summation_cf;

u_0 = ubar_cf.*(pop_start.^omega);
u_0 = u_0./geometric_mean(u_0);

error = 1000;

while error > 0.001

    kernel_utility_0 = u_0.*(lambda_cf.^(1/theta)); % The experience premium should appear as A_o^(1-betta) but does not matter, since it appears at the denominator too

    pi_0 = zeros(N,N,S);

    for j=1:N
        denom = sum(sum((kernel_utility_0./repmat(mu(j,:)',1,S)).^zeta));
        % Probability of going from j to any destination (:):
        pi_0(j,:,:) = ((kernel_utility_0./repmat(mu(j,:)',1,S)).^zeta)./denom;
    end
   
    L_y_cf = zeros(N,S);
    for j=1:N
        L_y_cf = L_y_cf + L_k_start(j)*squeeze(pi_0(j,:,:));
    end
    L_o_cf = L_y_start; L_k_cf = fertility_cf*sum(L_y_cf,2);
    pop_0 = L_k_cf + sum(L_y_cf,2) + sum(L_o_cf,2);
    
    u_1 = ubar_cf.*(pop_0.^omega);
    u_1 = u_1./geometric_mean(u_1);

    error = max(abs(u_0 - u_1));

    u_0 = 0.9*u_0 + 0.1*u_1;

end

pop_cf = pop_0;
pi_cf = pi_0;

end

